<!DOCTYPE html>
<html lang="en">
<head>
    <meta charset="UTF-8">
    <title>Fluid Simulation</title>

    <meta name="viewport" content="width=400">

    <link rel="stylesheet" type="text/css" href="dependencies/bootstrap.min.css">
    <link rel="stylesheet" type="text/css" href="dependencies/flat-ui.min.css">
    <link rel="stylesheet" type="text/css" href="main.css">

    <script>
      (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
      (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
      m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
      })(window,document,'script','https://www.google-analytics.com/analytics.js','ga');

      ga('create', 'UA-86531114-7', 'auto');
      ga('send', 'pageview');

    </script>

    <script id="2d-vertex-shader" type="x-shader/x-vertex">

        attribute vec2 a_position;

        void main() {
           gl_Position = vec4(a_position, 0, 1);
        }
    </script>

    <script id="boundaryShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_texture;
        uniform float u_scale;
        uniform vec2 u_textureSize;

        void main() {
            vec2 fragCoord = gl_FragCoord.xy;

            if (fragCoord.x < 1.0){
                gl_FragColor = u_scale*texture2D(u_texture, (fragCoord + vec2(1.0, 0.0))/u_textureSize);
                return;
            } else if (fragCoord.x >= u_textureSize.x-1.0){
                gl_FragColor = u_scale*texture2D(u_texture, (fragCoord + vec2(-1.0, 0.0))/u_textureSize);
                return;
            } else if (fragCoord.y < 1.0){
                gl_FragColor = u_scale*texture2D(u_texture, (fragCoord + vec2(0.0, 1.0))/u_textureSize);
                return;
            } else if (fragCoord.y >= u_textureSize.y-1.0){
                gl_FragColor = u_scale*texture2D(u_texture, (fragCoord + vec2(0.0, -1.0))/u_textureSize);
                return;
            }

            gl_FragColor = texture2D(u_texture, (fragCoord)/u_textureSize);
        }
    </script>

    <script id="2d-render-shader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_material;
        uniform vec2 u_textureSize;

        void main() {
            vec2 fragCoord = gl_FragCoord.xy;

            vec3 background = vec3(0.96, 0.87, 0.68);
            vec3 material1 = vec3(0.925, 0, 0.55);
            vec3 material2 = vec3(0.0, 0.70, 0.63);
            vec3 material3 = vec3(0.52, 0.81, 0.70);
            vec3 material4 = vec3(1.0, 0.7, 0.07);

            float val = texture2D(u_material, fragCoord/u_textureSize).x/2.0;
            if (val > 1.0) val = 1.0;
            if (val < 0.0) val = 0.0;

            float numColors = 3.0;

            vec3 color = vec3(0.0);
            if (val <= 1.0/numColors) {
                val *= numColors;
                color = background*(1.0-val) + material1*val;
            } else if (val <= 2.0/numColors) {
                val -= 1.0/numColors;
                val *= numColors;
                color = material1*(1.0-val) + material2*val;
            } else if (val <= 3.0/numColors) {
                val -= 2.0/numColors;
                val *= numColors;
                color = material2*(1.0-val) + material3*val;
            } else {
                val -= 3.0/numColors;
                val *= numColors;
                color = material3*(1.0-val) + material4*val;
            }

            gl_FragColor = vec4(color, 1);
        }
    </script>

    <script id="gradientSubtractionShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_velocity;
        uniform sampler2D u_pressure;

        uniform vec2 u_textureSize;

        uniform float u_const;

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;

            vec2 currentVelocity = texture2D(u_velocity, fragCoord/u_textureSize).xy;

            float n = texture2D(u_pressure, (fragCoord+vec2(0.0, 1.0))/u_textureSize).x;
            float s = texture2D(u_pressure, (fragCoord+vec2(0.0, -1.0))/u_textureSize).x;
            float e = texture2D(u_pressure, (fragCoord+vec2(1.0, 0.0))/u_textureSize).x;
            float w = texture2D(u_pressure, (fragCoord+vec2(-1.0, 0.0))/u_textureSize).x;

            gl_FragColor = vec4(currentVelocity-u_const*vec2(e-w, n-s), 0, 0);
        }
    </script>

    <script id="divergenceShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_velocity;

        uniform vec2 u_textureSize;

        uniform float u_const;

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;

            //finite difference formulation of divergence

            float n = texture2D(u_velocity, (fragCoord+vec2(0.0, 1.0))/u_textureSize).y;
            float s = texture2D(u_velocity, (fragCoord+vec2(0.0, -1.0))/u_textureSize).y;
            float e = texture2D(u_velocity, (fragCoord+vec2(1.0, 0.0))/u_textureSize).x;
            float w = texture2D(u_velocity, (fragCoord+vec2(-1.0, 0.0))/u_textureSize).x;

            float div = u_const*(e-w + n-s);
            gl_FragColor = vec4(div, 0, 0, 0);
        }
    </script>

    <script id="addMaterialShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_material;

        uniform vec2 u_textureSize;

        uniform vec2 u_mouseCoord;
        uniform float u_mouseLength;
        uniform float u_mouseEnable;

        uniform float u_reciprocalRadius;

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;

            float currentMaterial = texture2D(u_material, fragCoord/u_textureSize).x;

            if (u_mouseEnable == 1.0){
                vec2 pxDist = fragCoord - u_mouseCoord;
                currentMaterial += u_mouseLength*0.1*exp(-(pxDist.x*pxDist.x+pxDist.y*pxDist.y)*u_reciprocalRadius);
            }

            if (currentMaterial > 0.0) currentMaterial -= 0.002;//material disappears over time
            gl_FragColor = vec4(currentMaterial, 0, 0, 0);
        }
    </script>

    <script id="forceShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_velocity;

        uniform vec2 u_textureSize;

        uniform vec2 u_mouseCoord;
        uniform vec2 u_mouseDir;
        uniform float u_mouseEnable;

        uniform float u_reciprocalRadius;

        uniform float u_dt;

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;

            vec2 currentVelocity = texture2D(u_velocity, fragCoord/u_textureSize).xy;

            if (u_mouseEnable == 1.0){
                vec2 pxDist = fragCoord - u_mouseCoord;
                currentVelocity += u_mouseDir*u_dt*exp(-(pxDist.x*pxDist.x+pxDist.y*pxDist.y)*u_reciprocalRadius);
            }

            gl_FragColor = vec4(currentVelocity, 0, 0);
        }
    </script>

    <script id="jacobiShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_b;
        uniform sampler2D u_x;

        uniform vec2 u_textureSize;

        uniform float u_alpha;
        uniform float u_reciprocalBeta;

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;

            vec2 currentState = texture2D(u_b, fragCoord/u_textureSize).xy;

            //implicitly solve diffusion via jacobi iteration

            vec2 n = texture2D(u_x, (fragCoord+vec2(0.0, 1.0))/u_textureSize).xy;
            vec2 s = texture2D(u_x, (fragCoord+vec2(0.0, -1.0))/u_textureSize).xy;
            vec2 e = texture2D(u_x, (fragCoord+vec2(1.0, 0.0))/u_textureSize).xy;
            vec2 w = texture2D(u_x, (fragCoord+vec2(-1.0, 0.0))/u_textureSize).xy;

            vec2 nextState = (n + s + e + w + u_alpha * currentState) * u_reciprocalBeta;

            gl_FragColor = vec4(nextState, 0, 0);
        }
   </script>

    <script id="advectShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_velocity;
        uniform sampler2D u_material;

        uniform vec2 u_textureSize;
        uniform float u_scale;

        uniform float u_dt;

        vec2 bilinearInterp(vec2 pos, sampler2D texture, vec2 size){
            //bilinear interp between nearest cells

            vec2 pxCenter = vec2(0.5, 0.5);

            vec2 ceiled = ceil(pos);
            vec2 floored = floor(pos);

            vec2 n = texture2D(texture, (ceiled+pxCenter)/size).xy;//actually ne
            vec2 s = texture2D(texture, (floored+pxCenter)/size).xy;//actually sw
            if (ceiled.x != floored.x){
                vec2 se = texture2D(texture, (vec2(ceiled.x, floored.y)+pxCenter)/size).xy;
                vec2 nw = texture2D(texture, (vec2(floored.x, ceiled.y)+pxCenter)/size).xy;
                n = n*(pos.x-floored.x) + nw*(ceiled.x-pos.x);
                s = se*(pos.x-floored.x) + s*(ceiled.x-pos.x);
            }
            vec2 materialVal = n;
            if (ceiled.y != floored.y){
                materialVal = n*(pos.y-floored.y) + s*(ceiled.y-pos.y);
            }
            return materialVal;
        }

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;

            vec2 pxCenter = vec2(0.5, 0.5);

            vec2 currentVelocity;
            if (u_scale == 1.0) currentVelocity = 1.0/u_scale*texture2D(u_velocity, fragCoord/u_textureSize).xy;
            else {
                vec2 scaledCoord = (fragCoord-pxCenter)*u_scale;
                vec2 scaledSize = u_textureSize*u_scale;
                currentVelocity = 1.0/u_scale*bilinearInterp(vec2(1.0, 1.0) + scaledCoord/scaledSize*(scaledSize-vec2(0.5, 0.5)/u_scale), u_velocity, scaledSize);
            }

            //implicitly solve advection

            if (length(currentVelocity) == 0.0) {//no velocity
                gl_FragColor = vec4(texture2D(u_material, fragCoord/u_textureSize).xy, 0, 0);
                return;
            }

            vec2 pos = fragCoord - pxCenter - u_dt*currentVelocity;

            vec2 materialVal;
            //empty boundary
            if (pos.x < 0.0 || pos.x >= u_textureSize.x-1.0 || pos.y < 0.0 || pos.y >= u_textureSize.y-1.0) materialVal = vec2(0.0);
            else materialVal = bilinearInterp(pos, u_material, u_textureSize);

            gl_FragColor = vec4(materialVal, 0, 0);
        }
    </script>

    <script id="moveParticlesShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform sampler2D u_particles;
        uniform sampler2D u_velocity;

        uniform vec2 u_textureSize;
        uniform vec2 u_screenSize;
        uniform vec2 u_velocityTextureSize;

        uniform float u_dt;

        uniform float u_scale;

        vec2 bilinearInterp(vec2 pos, sampler2D texture, vec2 size){
            //bilinear interp between nearest cells

            vec2 pxCenter = vec2(0.5, 0.5);

            vec2 ceiled = ceil(pos);
            vec2 floored = floor(pos);

            vec2 n = texture2D(texture, (ceiled+pxCenter)/size).xy;//actually ne
            vec2 s = texture2D(texture, (floored+pxCenter)/size).xy;//actually sw
            if (ceiled.x != floored.x){
                vec2 se = texture2D(texture, (vec2(ceiled.x, floored.y)+pxCenter)/size).xy;
                vec2 nw = texture2D(texture, (vec2(floored.x, ceiled.y)+pxCenter)/size).xy;
                n = n*(pos.x-floored.x) + nw*(ceiled.x-pos.x);
                s = se*(pos.x-floored.x) + s*(ceiled.x-pos.x);
            }
            vec2 materialVal = n;
            if (ceiled.y != floored.y){
                materialVal = n*(pos.y-floored.y) + s*(ceiled.y-pos.y);
            }
            return materialVal;
        }

        void main() {

            vec2 fragCoord = gl_FragCoord.xy;
            vec2 particleCoord = texture2D(u_particles, fragCoord/u_textureSize).xy;

            vec2 currentVelocity = 1.0/u_scale*bilinearInterp(vec2(1.0, 1.0) + particleCoord*u_scale/u_velocityTextureSize*(u_velocityTextureSize-vec2(0.5, 0.5)/u_scale), u_velocity, u_velocityTextureSize);
            vec2 nextPosition = particleCoord+currentVelocity*u_dt;//explicitly solve advection

            if (nextPosition.x < 0.0) nextPosition.x = 0.0;
            else if (nextPosition.x >= u_screenSize.x-3.0) nextPosition.x = u_screenSize.x-3.0;
            if (nextPosition.y < 0.0) nextPosition.y = 0.0;
            else if (nextPosition.y >= u_screenSize.y-3.0) nextPosition.y = u_screenSize.y-3.0;

            gl_FragColor = vec4(nextPosition, 0, 0);
        }
    </script>

    <script id="packToBytesShader" type="x-shader/x-fragment">
        precision mediump float;

        uniform vec2 u_floatTextureDim;
        uniform sampler2D u_floatTexture;
        uniform float u_vectorLength;


        float shift_right (float v, float amt) {
            v = floor(v) + 0.5;
            return floor(v / exp2(amt));
        }
        float shift_left (float v, float amt) {
            return floor(v * exp2(amt) + 0.5);
        }
        float mask_last (float v, float bits) {
            return mod(v, shift_left(1.0, bits));
        }
        float extract_bits (float num, float from, float to) {
            from = floor(from + 0.5); to = floor(to + 0.5);
            return mask_last(shift_right(num, from), to - from);
        }
        vec4 encode_float (float val) {
            if (val == 0.0) return vec4(0, 0, 0, 0);
            float sign = val > 0.0 ? 0.0 : 1.0;
            val = abs(val);
            float exponent = floor(log2(val));
            float biased_exponent = exponent + 127.0;
            float fraction = ((val / exp2(exponent)) - 1.0) * 8388608.0;
            float t = biased_exponent / 2.0;
            float last_bit_of_biased_exponent = fract(t) * 2.0;
            float remaining_bits_of_biased_exponent = floor(t);
            float byte4 = extract_bits(fraction, 0.0, 8.0) / 255.0;
            float byte3 = extract_bits(fraction, 8.0, 16.0) / 255.0;
            float byte2 = (last_bit_of_biased_exponent * 128.0 + extract_bits(fraction, 16.0, 23.0)) / 255.0;
            float byte1 = (sign * 128.0 + remaining_bits_of_biased_exponent) / 255.0;
            return vec4(byte4, byte3, byte2, byte1);
        }

        void main(){
            vec2 fragCoord = gl_FragCoord.xy;
            float textureXcoord = floor((fragCoord.x - 0.5)/u_vectorLength+0.0001) + 0.5;
            vec4 data = texture2D(u_floatTexture, vec2(textureXcoord, fragCoord.y)/u_floatTextureDim);
            int textureIndex = int(floor(mod(fragCoord.x-0.5+0.0001, u_vectorLength)));
            if (textureIndex == 0) gl_FragColor = encode_float(data[0]);
            else if (textureIndex == 1) gl_FragColor = encode_float(data[1]);
            else if (textureIndex == 2) gl_FragColor = encode_float(data[2]);
            else if (textureIndex == 3) gl_FragColor = encode_float(data[3]);
        }
    </script>

    <script type="text/javascript" src="dependencies/jquery-3.1.0.min.js"></script>
    <script type="text/javascript" src="dependencies/flat-ui.min.js"></script>
    <script type="text/javascript" src="dependencies/three.min.js"></script>

    <script type="text/javascript" src="js/threeView.js"></script>
    <script type="text/javascript" src="js/GLBoilerplate.js"></script>
    <script type="text/javascript" src="js/GPUMath.js"></script>
    <script type="text/javascript" src="js/main.js"></script>
</head>
<body>

<canvas id="glcanvas"></canvas>
<div id="threeContainer"></div>

<a href="#" id="about">?</a>

<div class="modal fade" id="aboutModal" tabindex="-1" role="dialog" aria-labelledby="basicModal" aria-hidden="true">
    <div class="modal-dialog modal-lg">
        <div class="modal-content">
            <div class="modal-body">
                <b>Fluid Simulation Shader</b><br/><br/>
                This simulation solves the <a href="https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations" target="_blank">Navier-Stokes equations</a> for incompressible fluids in a GPU fragment shader.
                I implemented <a href="https://en.wikipedia.org/wiki/No-slip_condition" target="_blank">no-slip boundary conditions</a> at the borders to keep the fluid contained within the bounds of the screen.
                To increase performance, I solved for the velocity vector field of the fluid at a lower resolution than I used to compute the visualization of fluid flow; I used bilinear interpolation to smooth out artifacts caused by this speedup.
                I've also added 160,000 <a href="https://en.wikipedia.org/wiki/Lagrangian_particle_tracking" target="_blank">Lagrangian particles</a> on top of the simulation -
                these particles are rendered using <a href="https://threejs.org/" target="_blank">threejs</a>, but their positions are computed on the GPU.
                <br/><br/>
                <b>Instructions:</b> Click and drag to apply a force to the fluid.  Over time, the colored material in the fluid will dissipate:
                <img style="width:100%" src="img2.jpg"/>
                <br/><br/>
                To learn more about the math involved, check out the following sources:<br/>
                <a href="https://pdfs.semanticscholar.org/84b8/c7b7eecf90ebd9d54a51544ca0f8ff93c137.pdf" target="_blank">Real-time ink simulation using a grid-particle method</a> - mixing Eulerian and Lagrangian techniques for fluids.<br/>
                <a href="http://developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch38.html" target="_blank">Fast Fluid Dynamics Simulation on the GPU</a> - a very well written tutorial about programming the Navier-Stokes equations on a GPU.
                Though not WebGL specific, it was still very useful.<br/>
                <a href="http://jamie-wong.com/2016/08/05/webgl-fluid-simulation/" target="_blank">Fluid Simulation (with WebGL demo)</a> - this article has some nice, interactive graphics that helped me debug my code.<br/>
                <a href="http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/ns.pdf" target="_blank">Stable Fluids</a> - a paper about stable numerical methods for evaluating Navier-Stokes on a discrete grid.<br/>
                <br/>
                By <a href="http://www.amandaghassaei.com/" target="_blank">Amanda Ghassaei</a>, code on <a href="https://github.com/amandaghassaei/FluidSimulation" target="_blank">Github</a>.
                <br/><br/>
            </div>
        </div>
    </div>
</div>

<div class="modal fade" id="noSupportModal" tabindex="-1" role="dialog" aria-labelledby="basicModal" aria-hidden="true">
    <div class="modal-dialog modal-lg">
        <div class="modal-content">
            <div class="modal-body">
                This app uses GPU functions that are not supported by this device/browser, please try again on desktop with Chrome/Firefox.<br/><br/>
            </div>
        </div>
    </div>
</div>

</body>
</html>